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ABSTRACT 

We present the results of high-resolution, three-dimensional (3D) hydrody- 
namic simulations of the dynamics and formation of coherent, long-lived vortices 
in stably-stratified protoplanetary disks. Tall, columnar vortices that extend ver- 
tically through many scale heights in the disk are unstable to small perturbations; 
such vortices cannot maintain vertical alignment over more than a couple scale 
heights and are ripped apart by the Keplerian shear. Short, finite-height vortices 
that extend only one scale height above and below the midplane are also unsta- 
ble, but for a different reason: we have isolated an antisymmetric (with respect 
to the midplane) eigenmode that grows with an e-folding time of only a few or- 
bital periods; the nonlinear evolution of this instability leads to the destruction 
of the vortex. Serendipitously, we observe the formation of 3D vortices that are 
centered not in the midplane, but at one to three scale heights above and below. 
Breaking internal gravity waves create vorticity; anticyclonic regions of vorticity 
roll-up and coalesce into new vortices, whereas cyclonic regions shear into thin 
azimuthal bands. Unlike the midplane-centered vortices that were placed ad hoc 



1 NSF Astronomy & Astrophysics Postdoctoral Fellow 



- 2 - 



in the disk and turned out to be linearly unstable, the off-midplane vortices form 
naturally out of perturbations in the disk, and are stable and robust for many 
hundreds of orbits. 

Subject headings: accretion, accretion disks — hydrodynamics — planetary sys- 
tems: protoplanetary disks 

1. INTRODUCTION 

Like the atmosphere of Jupiter, protoplanetary disks are characterized by rapid rotation 
and intense shear, inspiring proposals that disks may also be populated with long-lived, ro- 
bust storms analogous to the Great Red Spot (Barge & Sommeria 1995; Adams & Watkins 
1995). Such vortices may play key roles in the formation of stars and planets: (1) in cool, 
neutral protoplanetary disks, vortices might transport angular momentum radially outward 
so that mass can continue to accrete onto the growing protostar; (2) models of protoplanet 
migration (that account for "Hot Jupiters," gas giant planets that closely orbit their parent 
stars with periods of only a few days) are highly sensitive to the effective "viscosity" or 
turbulent transport within disks; transport mediated by coherent vortices rather than tur- 
bulent eddies could lead to qualitative and quantitative changes to migration rates; and (3) 
vortices rapidly sweep-up and concentrate dust particles, which may help in the formation 
of kilometer-size planetesimals, the "building blocks" of planets, either by increasing the 
efficiency of binary agglomeration or by triggering a local gravitational instability. 

1.1. The angular momentum problem in accretion disks 

One of the most vexing problems in astrophysics is how mass and angular momentum 
are transported in accretion disks. Because disks are differentially rotating, one might think 
that viscous torquing could transport angular momentum outward and mass inward onto the 
central object. However, if the source of viscosity is molecular collisions, then the time-scale 
for such viscous transport would exceed the age of the universe by many orders of magnitude. 
Shakura & Sunyaev (1973) proposed that turbulence might enhance the transport, and so 
introduced a turbulent or "eddy" viscosity: u t = aHc s , where H is the pressure scale height 
of the disk, c s is the sound speed, and a is a dimensionless number (less than unity because 
turbulent eddies would most likely be subsonic and no larger than the scale height of the 
disk). The source of this turbulence (e.g., shear instabilities, convection, finite-amplitude 
perturbations, magnetic fields), to say the least, has been highly controversial. See Balbus 
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& Hawley (1998) and Stone et al. (2000) for thorough reviews. 

Balbus & Hawley (1991) (also Hawley & Balbus 1991) applied the magneto-rotational 
instability (MRI) of Velikhov (1959) and Chandrasekhar (1960, 1961), and demonstrated 
that weak magnetic fields can destabilize the Keplerian shear in an accretion disk, leading to 
the creation of turbulence and the outward transport of angular momentum. However, the 
relatively cool and nearly neutral accretion disks around protostars most likely lack sufficient 
coupling between matter and magnetic fields (Blaes & Balbus 1994), except perhaps in thin 
surface layers that have been ionized by cosmic rays or protostellar x-rays (Gammie 1996). 

Thus, there still remains strong interest in finding a purely hydrodynamic means of mass 
and angular momentum transport. Bracco et al. (1998), using two-dimensional, incompress- 
ible fluid dynamics, have shown that long-lived, coherent anticyclonic vortices emerge out 
of a Keplerian shear flow that is seeded with small-scale vorticity perturbations; they ob- 
served smaller vortices merging to form larger vortices, demonstrating the "inverse cascade" 
of energy from small to large scales that is a hallmark of two-dimensional turbulent flows 
(Provenzale 1999). Godon & Livio (1999, 2000) confirmed this result with two-dimensional, 
compressible, barotropic simulations. Lovelace et al. (1999) have found a linear instability of 
nonaxisymmetric Rossby waves in thin, non- magnetized Keplerian disks when there is a local 
extremum in an entropy-modified potential vorticity. Li et al. (2000, 2001) showed that such 
Rossby waves break and coalesce to form vortices which radially migrate and transport mass 
through the disk. Klahr & Bodenheimer (2003) have demonstrated that a globally unstable 
radial entropy gradient in a protoplanetary disk generates Rossby waves, vortices and turbu- 
lence that may result in local outward transport of angular momentum. To date, there are 
still unresolved issues regarding realistic formation mechanisms and long-term maintenance 
of vortices in disks. 



1.2. The role of turbulence and vortices in planet formation 

Whether protoplanetary disks are filled with turbulent eddies or long-lived coherent 
vortices (or both) is critical for not only accretion, but also for two key processes in planet 
formation: protoplanet migration and planetesimal formation. Extrasolar planet surveys 
reveal a class of gas giant planets that closely orbit their parent stars with periods of only 
a few days (Mayor & Queloz 1995; Marcy & Butler 1998; Marcy et al. 2000). These "Hot 
Jupiters" most likely did not form in situ, but instead formed farther out in their disks where 
it was cooler, and then migrated inward to their present locations (Lin et al. 1996). The exact 
mechanism for this migration depends on the size of the protoplanet: a small protoplanet 
(< 10M e ) raises tides in the disk which exert torques back on the protoplanet, causing it to 
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migrate inward (Type I migration), whereas a larger protoplanet opens up a gap in the disk, 
and both gap and protoplanet migrate inward together on a slow "viscous" timescale (Type 
II migration) (Ward 1997). In numerical simulations of protoplanet migration, the viscous 
evolution is usually modeled with an eddy viscosity. Although such a crude turbulence model 
is useful for estimating global disk properties, it is not at all clear that it is appropriate to 
apply it to complex disk-planet interactions. Does a turbulent disk really behave exactly the 
same way as a laminar disk with just a larger viscosity? Qualitatively and quantitatively, 
transport and mixing due to long-lived coherent vortices are not necessarily well modeled by 
an enhanced viscosity (Roller et al. 2003). 

Another aspect of planet formation that is extremely sensitive to turbulence is the forma- 
tion of planetesimals, the kilometer-size "building blocks" of planets. In a quiescent, laminar 
protoplanetary disk, planetesimals may form directly from the gravitational clumping of a 
thin, dense dust sublayer that has settled into the midplane (Safranov 1960; Goldreich & 
Ward 1973). However, turbulence (either from MRI, or generated by Kelvin-Helmholtz-like 
instabilities of the vertical shear between the dust sub-layer in the midplane, which orbits 
at the Keplerian velocity, and the gas-rich layers immediately above and below it, which 
orbit at a slightly slower rate because of the outward gas pressure force) might "kick up" the 
dust and prevent the sub-layer from reaching the critical density necessary for gravitational 
instability (Cuzzi et al. 1993; Weidenschilling 1995; Champney et al. 1995). Further growth 
of grains must instead continue via binary agglomeration (Weidenschilling & Cuzzi 1993). 

Vortices may aid in the formation of planetesimals in either scenario. Barge & Sommeria 
(1995), Tanga et al. (1996), Chavanis (2000), and de la Fuente Marcos & Barge (2001) 
tracked the motion of individual particles in analytic models of two-dimensional anticyclones 
embedded in Keplerian shear. Klahr & Henning (1997) considered simple models of vortices 
that corresponded to convective cells. Johansen et al. (2004) used two-fluid simulations to 
study the interactions of the gas and dust components inside vortices. However, because 
of the relatively low resolution and high numerical dissipation in their simulations, their 
vortices were short-lived and decayed away in roughly one orbit. In all these models, gas 
drag caused grains to spiral into the vortex centers. Density enhancements inside vortices 
might increase the efficiency of binary agglomeration, or might trigger local gravitational 
clumping within the vortex center as opposed to the global gravitational clumping of the 
whole dust sub-layer envisioned in the original Safranov (1960) and Goldreich & Ward (1973) 
theory. 
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1.3. Outline 

This work is the first in a series on high-resolution numerical simulations of the dynamics 
and formation of long-lived, coherent vortices in protoplanetary disks, and on what role they 
may play in star and planet formation. In §2, we will overview dimensional analysis and 
key timescales for vortices in protoplanetary disks. We will discuss the deficiencies of 2D 
analyses and explain the need for 3D simulations. In §3, we will present the hydrodynamical 
equations for the base disk and the vortex flow, as well as briefly describe the numerical 
method. Results of numerical simulations will be presented in §4, including a discussion on 
angular momentum transport. Conclusions and avenues for future work appear in §5. 



2. A PRIMER FOR VORTEX DYNAMICS IN DISKS 

2.1. Dimensional analysis and key timescales 

Protoplanetary disks are cool in the sense that the gas sound speed c s is much slower 
than the Keplerian orbital velocity Vk{t) = rf2^;(r) = a/ GM+/r, where G is the gravita- 
tional constant, M* is the protostellar mass, and r is the cylindrical radius to the protostar. 
Hydrostatic balance implies that the time it takes sound waves to traverse the thickness of 
the disk is of order the orbital period. Thus, cool disks are geometrically thin (see Frank 
et al. 1985): 

c s ~ Q K H, (2-la) 
S = c s /V K ~ H/r<l, (2-lb) 

where H is the vertical pressure scale height. The radial component of the protostellar 
gravity nearly balances the centrifugal force, but because of the relatively weak outward 
radial pressure force, the gas orbits at slightly slower than the Keplerian velocity: 

Q(r,z) = Q K (r)[l- V (r,z)], (2-2) 

where rj ~ 0(5 2 ). The horizontal shear rate of this base flow is: 

_ dQ K 3 
°k = r— = --n K . (2-3) 

Keplerian shear is anticyclonic (as indicated by the negative sign), and is comparable in 
magnitude to the rotation rate itself. As we will see, this has important consequences for 
the properties of vortices in protoplanetary disks. 
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Shear will stretch and tear a vortex apart unless: (1) the vortex rotates in the same 
sense as the ambient shear, and (2) the vortex is at least as strong as the shear (Moore & 
Saffman 1971; Kida 1981; Marcus 1990, 1993). The first condition implies that all long-lived 
vortices in a protoplanetary disk will be anticyclones; that is, the vortices will rotate in 
a sense opposite of that of the overall rotation of the disk. Let v± be the characteristic 
horizontal speed of gas around a vortex with aspect ratio x = A^,/A r > 1, where A r (along 
the radial direction in the disk) is the minor axis and (along the azimuthal direction 
in the disk) is the major axis. The characteristic ^-component of vorticity 1 in the core of 
the vortex is Cj z ~ v±/A r . The second condition above implies that the horizontal speed 
of gas around the vortex must be of order the differential velocity across the vortex due 
to the ambient shear, or equivalently, that the characteristic vorticity is of order the shear 
rate: uj z ~ v±/A r ~ <r^. Moore & Saffman (1971) and Kida (1981) analytically determined 
the exact steady-state solution for a 2D elliptical vortex of uniform vorticity embedded in a 
uniform shear and showed that: 

5i = i (2-4) 

for a vortex that rotates in the same sense as the shear. A vortex that is weak relative 
to the shear will be stretched out and have a large aspect ratio, whereas a vortex that is 
comparable in strength to the shear will be more compact and have an aspect ratio closer 
to unity. Because the Keplerian shear is comparable in magnitude to the Keplerian rotation 
rate, the rotational period of gas around the core of the vortex, r vor = 4ir/uj z , is of order 
the orbital period of the vortex around the protostar, r or & = 2n/VL. The Rossby number is 
defined to be the ratio of these two timescales: 

_ Cb g r orb 3 x + 1 1 (n -x 

Ro = — ~ ~ -— , (2-5) 

2Q r vor Ax-lx 

which is order unity for compact vortices. In contrast, the Great Red Spot (GRS) on Jupiter 
rotates with a period of roughly six days, whereas the planet rotates in just under 10 hours, 
implying a Rossby number Ro xs 0.18 (Marcus 1993). 2 

As we are interested in long-lived vortices, we require that the gas velocity be subsonic; 
otherwise, sound waves and shocks would rapidly dissipate the kinetic energy of the vortex. 
Equations (2-la) and (2-5) imply: 

e = ^ ~ ±Ro, (2-6) 
C„ ti 



1 Recall that vorticity is defined u = V x v, and is two times the local angular velocity of the fluid. 

2 For a rotating planet, Ro ~ (r or (,/sin X)/T vor , where A is the latitude. Only the component of the 
rotation normal to the surface of the planet enters into the horizontal component of the Coriolis force. 
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where e is defined to be the horizontal Mach number. Thus, the horizontal extent of subsonic, 
compact vortices in a Keplerian shear cannot be much greater than the scale height of the disk. 
Such vortices are 3D in nature. In contrast, the GRS has considerably more "pancake-like" 
dimensions: 26,000 km x 13,000 km x 40 km (Marcus 1993). 

Another important timescale is the Brunt-Vaisala period: tbv = Zk/ubv, which is the 
period of oscillations in a convectively stable atmosphere (e.g., internal gravity waves). The 
Brunt-Vaisala frequency is (see Pedlosky 1979): oobv = a/ (g/C p )ds/dz, where g is the local 
gravitational acceleration, C p is the specific heat at constant pressure, and s is the mean 
entropy. If the protoplanetary disk is vertically isothermal, then lvbv ~ Qk\z\/H. That is, 
except in the immediate vicinity of the midplane, the Brunt-Vaisala period is also of order 
the orbital period in the disk. The internal Froude number is defined as the ratio of the 
Brunt-Vaisala period to the rotational period of the gas in a vortex: 



When the Froude number is much less than unity, a system is considered to be strongly 
stratified; large-scale vertical motions are suppressed and the different layers of the fluid are 
only weakly coupled. 

Table 1 summarizes the values of the key timescales discussed so far. For comparison, 
values for the GRS are presented as well. On Jupiter, these timescales are well ordered: 
r vor ^> T or b ^> tbv] whereas in a protoplanetary disk: r vor ~ T or b ~ tbv- The near-equality 
of timescales in a protoplanetary disk simply follows from the fact that at any given location 
in the disk, the protostellar gravity sets the orbital velocity, shear, and strength of the 
stratification. Because the Brunt-Vaisala frequency for the GRS is so much faster than the 
other timescales, the vertical and horizontal dynamics can be decoupled; the GRS is well- 
modeled with two-dimensional fluid dynamics. No such clear separation of timescales occurs 
for a protoplanetary disk, and the horizontal and vertical dynamics cannot be decoupled. 
In general, when the Rossby and Froude numbers are small (e.g., rotation and stratification 
dominate over shear and nonlinear advection), then the horizontal and vertical motions can 
usually be decoupled and two-dimensional analyses can safely be applied. 

In many of the common 2D models (e.g., "shallow- water" , quasigeostrophic, or 2D- 
barotropic), vorticity (or a potential vorticity) is an advectively conserved quantity (because 
baroclinicity and vortex tilting are neglected). In 2D, only the ^-component of vorticity is 
nonzero, and vortex lines can only be oriented vertically. The rotation axes of long-lived 
coherent vortices as well as transient turbulent eddies are constrained to be aligned parallel 
(or anti-parallel) to the vertical axis. Furthermore, shear will stretch out and destroy those 
vortices that do not have the same sense of rotation as the shear. Thus, a 2D shear flow 




(2-7) 
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will be populated with vortices that are all perfectly aligned vertically and with the same 
sense of rotation. It has been demonstrated in both laboratory experiments and numerical 
simulations that in such restricted flows, small regions of vorticity merge to form larger 
vortices — this phenomenon is called an "inverse cascade" of energy from small to large 
scales (Kraichnan 1967; Lesieur 1997; Paret & Tabeling 1998; Baroud et al. 2003) Out of an 
apparently chaotic flow filled with small-scale, transient eddies, large-scale coherent features 
can emerge. 

In 3D, the dynamics of turbulence are very different. Vortex lines can be oriented in any 
direction, and they can twist and bend. Vortices and eddies tilt and stretch their neighbors, 
disrupting them and eventually destroying them as they give up their energy to smaller 
eddies. In fact, this is the foundation for the Richardson and Kolmogorov model of 3D, 
isotropic turbulence: turbulent kinetic energy is transfered (via nonlinear interactions) from 
large to small to smaller eddies on down (i.e., "forward cascade") until it is destroyed by 
viscous dissipation. 

An open question (and a very active area of basic research in fluid dynamics) is what 
happens when rotation and stratification are important, but not overwhelmingly dominant 
(e.g., when the Rossby and/or Froude numbers are of order unity)? How much rotation and 
stratification are necessary for a real, 3D flow to exhibit 2D characteristics, such as inverse 
cascades? Bracco et al. (1998) and Godon & Livio (1999, 2000) showed that small vortices 
merged to form larger vortices in 2D simulations, but these results have yet to be confirmed 
with 3D simulations. 



2.2. A model of a 3D vortex 

One of the best ways to visualize a vortex in three dimensions is to graph vortex lines, 
curves that are everywhere tangent to the vorticity vector field. Vorticity is, by definition, 
divergence-free; thus, vortex lines (like magnetic field lines) cannot have beginnings or end- 
ings within the fluid, but must extend to the boundaries, or off to infinity, or form closed 
loops. The simplest example of a 3D vortex is an infinite column of rotating fluid in which 
the fluid velocity is independent of height. The core of uniform vorticity is threaded by 
infinitely-long parallel vortex lines (see Figure la). One might imagine chopping off the ends 
of an infinite column to create a finite-height cylinder of rotating fluid. However, vortex lines 
cannot have loose ends in the fluid, and instead must wrap-around and form closed loops 
(see Figure lb). Note that in the core of such a vortex, the vortex lines will be oriented in 
one direction, whereas in a halo surrounding the core, the vortex lines will be oriented in the 
opposite direction. 
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Timescale 


GRS 


PPD 


T vor = 4tt/uj z 


8 d. 


~ 1 y. 


r orh = 2n/Vt 


10 h. 


~ 1 y. 


TBV = ^/^BV 


6 m. 


~ 1 y. 


Ro = 1 7~yor 


0.18 


0(1) 


Fr = tbv/tvot 


5 x 10~ 4 


0(1) 


Ri = 1/Fr 2 


4 x 10 6 


0(1) 



Table 1: Comparison of important dynamical timescales for the Great Red Spot and a vortex 
in a protoplanetary disk. Here, we use "year" in the more general sense to refer to the 
orbital period in the disk. r vor is the orbital period of gas around the core of the vortex, 
r orb is the rotational period of the system, and r BV is the Brunt- Vaisala period, or the 
timescale for vertical thermal oscillations in a stably stratified atmosphere. Note that for 
the GRS, r vor ^> r or b ^> t B v, whereas for the PPD, r vor ~ T or b as tbv- The lower half of the 
table expresses the same data in terms of common nondimensional parameters from fluid 
dynamics: the Rossby number, the internal Froude number, and the Richardson number. 
For a planetary atmosphere, only the component of the rotation vector normal to the surface 
is relevant, so the Rossby number is larger by a factor 1/ sin A, where A is the latitude on 
the planet. 



Variable Scaling 

v±/c s e ~ (k r /H)Ro 
P/P ~ PIP ~ T/f e 2 /Ro 

v z /c s < e 2 



Table 2: How velocity and thermodynamic variables scale with Mach number e and Rossby 
number Ro. 



Fig. 1. — Vortex lines in an infinite columnar vortex and a finite- height cylindrical vortex. 
Downward oriented arrows denote anticyclonic vorticity, whereas upward-oriented arrows 
indicate cyclonic vorticity. 
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Fig. 2. — Balance of forces in an anticyclone with Rossby number less than unity. 
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We now turn to the balance of forces in a 3D vortex. In any horizontal plane, the 
centrifugal force always points radially outward from the vortex center, whereas the Coriolis 
force points outward for cyclones and inward for anticyclones. When the Rossby number is 
much greater than unity, the Coriolis force is negligible. The outward centrifugal force must 
be balanced by an inward pressure force; such vortices must have low-pressure cores. When 
the Rossby number is less than unity, the centrifugal force is small and the Coriolis force 
is balanced by the pressure force. Cyclones must have low-pressure cores and anticyclones 
must have high-pressure cores. 

In the vertical direction, the only force that can balance the pressure force is buoyancy. 
Figure 2a shows the balance of forces in a low Rossby number anticyclone that is vertically 
centered on the midplane, whereas figure 2b shows the same for an anticyclone located 
completely above the midplane. In the first case, the high-pressure anticyclone must have 
cool, dense lids to provide a buoyant force toward the midplane that balances the pressure 
force away from the midplane. In the second case, the top lid (i.e., the one farthest from 
the midplane) must be cool and dense, whereas the bottom lid (i.e., the one closest to the 
midplane) must be warm and light. These thermal lids will be weakened or destroyed either 
by radiation of heat or turbulent erosion. Thus, a vertical flow must develop within the 
vortex to maintain the lids. For a high-pressure anticyclone centered on the midplane, gas 
will rise away from the midplane, adiabatically expanding and cooling. At the lids, the gas 
diverges outward and recirculates back toward the midplane along the edges of the vortex. 
For an anticyclone that is completely above the midplane, gas will rise (sink) toward the top 
(bottom) lid, where it will expand and cool (compress and heat-up). 

We can use the arguments in this section to derive simple scaling relationships for the 
vortex flow variables, which will then guide us in making further approximations to the 
fluid equations. We decompose each flow variable into a time-independent, axisymmetric 
component describing the base disk flow, which we will denote with overbars, and a time- 
dependent component describing the vortex flow, which we will denote with tildes (e.g., 
p = p + p,p = p + p). As discussed in the previous subsection, we assume that the horizontal 
vortex flow is subsonic: v±/c s ~ e < 1. If the Rossby number is of order unity or less, then the 
horizontal pressure force must be of the same order as the Coriolis force (geostrophic balance): 
p/pA r ~ 2flv_ L , where p is the pressure perturbation associated with the vortex, and p is the 
mean density. Using (2- la) and (2-6), one can show that the pressure perturbation in the 
vortex must scale as: p/p ~ e 2 /Ro. We also assume the vertical buoyancy force is of order 
the vertical pressure gradient: g z T/T ~ p/pH, which leads to the temperature fluctuation 
having the same scaling as the pressure fluctuation: T/T ~ p/p. Because of the ideal 
gas law, this implies that density fluctuations also have the same scaling. The scaling for 
the vertical velocity is set by the maintenance of temperature perturbations in the thermal 
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lids. In the temperature equation, temperature fluctuations are created via pressure-volume 
work: v z T(ds/dz)/C v , where v z is the vertical velocity, s is the mean entropy profile, and 
C v is the heat capacity at constant volume. The destruction process is either radiative 
heating/cooling, T/r coo i, where t coo \ is a characteristic cooling time; or turbulent advection, 
T I Teddy ~ v±T/A r , where r e ddy is a characteristic eddy diffusion time. If the vortex flow is 
smooth and laminar, then we balance pressure-volume work with thermal cooling to obtain 
the scaling: v z ~ (e 2 /Ro)H/r coo i. Thus, if the vortex is laminar and the cooling time is 
long, the vertical velocity can be quite small. On the other hand, if the vortex flow is highly 
turbulent, we balance pressure- volume work and the turbulent advection of heat to obtain 
the scaling: v z ~ e 2 . These scaling relationships are summarized in Table 2. 



3. FLUID EQUATIONS AND NUMERICAL METHOD 

3.1. Cartesian and anelastic approximations 

Motivated by the scalings discussed in the previous section, we focus on subsonic, 
compact vortices that have horizontal extent A r comparable to the vertical scale height 
H, which is much smaller than the distance r to the protostar: A r /H ~ e/Ro < 1 and 
A r /r ~ e5/Ro <C 1. This allows us to make two key approximations to the hydrodynamics. 

The first simplification we make is the Cartesian approximation (Goldreich & Lynden- 
Bell 1965): we simulate the hydrodynamics only within a small patch of the disk (Ar <C 
To, A0 <C 27r) that co-rotates with the gas at some fiducial radius r . We map this patch of 
the disk onto a Cartesian grid: r — ro — > x, ro(0 — <f>o) — > y, z — > z, v r — > v x , v<f, — > v y , and 
v z —> Vz- The background shear and mean thermodynamic variables have radial gradients 
that vary on the length scale r which is much larger than the characteristic size of a subsonic 
vortex. For example, let q represent any background disk variable; then the variation Sq 
over the size of a vortex is: 5q/q ~ (<91ng/<9r)A r ~ A r /r I. This allows us to neglect 
radial gradients of the mean thermodynamic variables and to linearize the Keplerian shear 
flow. The time-independent, axisymmetric base Keplerian flow in the rotating frame, which 
we denote with overbars, is then: 



v = -|fi xy, (3-la) 

f = T , (3-lb) 

p = p ex P (-z 2 /2H*), (3-lc) 

p = p exp(-z 2 /2H 2 ), (3-ld) 

s = Kz 2 /2H%, (3-le) 
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where T, p, p, s are the mean temperature, density, pressure, and entropy, respectively, 
which depend only on the vertical coordinate z. We have defined: Q = Q K (r ), p = p lZT , 
Hq = TZTq/Qq, and TZ is the gas constant. The errors in making this Cartesian approximation 
are of order 0(eS, S 2 ). Note that neglecting radial gradients of the background disk properties 
effectively filters out Rossby waves (which require a gradient in the background vorticity) 
and large-scale baroclinic instabilities. Our rationale for this is that if there were indeed 
large-scale baroclinic instabilities, they would produce large, supersonic vortices which would 
rapidly decay from radiation of acoustic waves and shocks. 

The second simplification we make is the anelastic approximation: all variables are 
expanded in powers of the Mach number (e.g., v = e(v + v) + ... ; p = p + e 2 p + ... ; 
p = p + e 2 p + ... ). These expansions are then substituted into the Euler equations, and 
terms are grouped by like powers of e; we keep only the first-order corrections (Barranco 
et al. 2000). The anelastic approximation has been used extensively in the study of deep, 
subsonic convection in planetary atmospheres (Ogura & Phillips 1962; Gough 1969) and stars 
(Gilman & Glatzmaier 1981; Glatzmaier & Gilman 1981a,b). One of the consequences of this 
approximation is that the total density is replaced by the time-independent mean density in 
the mass continuity equation, which has the effect of filtering high-frequency acoustic waves 
and shocks, but allowing slower wave phenomena such as internal gravity waves. The errors 
in making the anelastic approximation are of order 0(e 2 ). The anelastic Euler equations in 
the co-rotating frame are: 

p = pUf + p-RT, (3-2a) 
= V-(pv), (3-2b) 

= -2Q z x v + 3Q 2 xit - -Vp - ^Q 2 zz, (3-2c) 
at p p 

^ = -(7-l)T(V.v)-(T-f)/r cooi (^), (3-2d) 

where the advective derivative is defined as d/dt = d/dt + v • V, and where 7 = Cp/Cy is 
the ratio of specific heats. Equation (3- 2b) is the crux of the anelastic approximation: the 
time-derivative has dropped out of the mass continuity equation, which effectively filters out 
all acoustic phenomena. The terms on the right hand side of the momentum equation (3-2c) 
are: the Coriolis force, the tidal force (the remainder after balancing the radial component 
of the protostellar gravity and the centrifugal non-inertial force), the pressure force, and the 
buoyancy force. In equation (3-2d), the first term on the right hand side is pressure- volume 
work (which is reversible) and the second term represents a simplified model of cooling 
or thermal relaxation. In the limit where the cooling time t coo \ is infinitely long, one can 
show that the temperature equation is equivalent to entropy being advectively conserved: 
ds/dt = 0. 
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3.2. Spectral expansions and boundary conditions 

In this and the following subsection, we briefly describe the numerical method; a more 
detailed presentation can be found in Barranco & Marcus (2004). We solve the Euler equa- 
tions (3-2) with a spectral method; that is, each variable is represented as a finite sum of basis 
functions multiplied by spectral coefficients (Gottlieb & Orszag 1977; Marcus 1986; Canuto 
et al. 1988; Boyd 1989). The choice of basis functions for each direction is guided by the cor- 
responding boundary conditions. The equations are autonomous in the azimuthal coordinate 
y, so it is reasonable to assume periodic boundary conditions in this direction. However, the 
equations explicitly depend on the radial coordinate x because of the linear background shear. 
We adopt "shearing box" boundary conditions: q(x + L x , y — (3/2)Q L x t, z, t) = q(x, y, z, t), 
where q represents any of v, p, T, etc. In practice, we rewrite the equations (3-2) in 
terms of quasi-Lagrangian or shearing coordinates that advect with the background shear 
(Goldreich & Lynden-Bell 1965; Marcus & Press 1977; Rogallo 1981): t' = t, x' = x, 
y' = y + (3/2)Q xt, and z' = z. In these new coordinates, the radial boundary conditions 
become: q(x' + L x ,y',z',t') = q(x' , y' , z' , t') . That is, shearing box boundary conditions 
are equivalent to periodic boundary conditions in the shearing coordinates. Physically, this 
means that the periodic images at different radii are not fixed, but advect with the back- 
ground shear. 

In the shearing coordinates, the equations are autonomous in both x' and y' (although 
they now explicitly depend on t'), making a Fourier basis the natural choice for the spectral 
expansions in the horizontal directions: 

q(x', y', z', f ) = ftMe^^), (3-3) 

k 

where q is any variable of interest, {qk(t')} is the set of spectral coefficients, and k = 
{k' x , k' y , n} is the set of wavenumbers. We have implemented the simulations with two differ- 
ent sets of basis functions for the spectral expansions in the vertical direction, corresponding 
to two different sets of boundary conditions: 

(i) For the truncated domain — L z < z < L z , we use Chebyshev polynomials: (f) n (z) = 
T n (z/L z ) = cos(n£), where £ = cos^ 1 (z/L z ). We apply the condition that the vertical 
velocity vanish at the top and bottom boundaries: v z (x,y, z = ±L z ,t) = 0. 

(ii) For the infinite domain — oo < z < oo, we use rational Chebyshev functions: (f> n (z) = 
cos(n£) (for v x , v y , u z , and all the thermodynamic variables) or (f> n (z) = sin(n£) (for v z , Cj x 
and 0J y ), where £ = cot" 1 {z / L z ) . In this context, L z is no longer the physical size of the 
box, but is a mapping parameter; exactly one half of the grid points are within \z\ < L z , 
whereas the other half are widely spaced in the region L z < \z\ < oo. No explicit boundary 
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conditions on the vertical velocity are necessary when we solve the equations on the infinite 
domain because the (f> n (z) = sin(n£) basis functions individually decay to zero at large z. 

3.3. Brief description of the numerical method 

The equations are integrated forward in time with a fractional step method: the nonlin- 
ear advection terms are integrated with an explicit second-order Adams-Bashforth method, 
and the pressure step is computed with a semi-implicit second-order Crank-Nicholson method. 
The time integration scheme is overall globally second-order accurate. Unlike finite-difference 
methods, spectral methods have no inherent grid dissipation; energy cascades to smaller and 
smaller size scales via the nonlinear interactions, where it can "pile-up" and potentially de- 
grade the convergence of the spectral expansions. We employ a V 8 hyperviscosity or low-pass 
filter every timestep to damp the energy at the highest wavenumbers. 

Different horizontal Fourier modes interact only through the nonlinear advective terms; 
once these terms are computed, the horizontal Fourier modes can be decoupled. This mo- 
tivated us to parallelize the code: each processor computes on a different block of data in 
horizontal Fourier wavenumber space. Parallelization is implemented with Message Passing 
Interface (MPI) on the IBM Blue Horizon and IBM Datastar at the San Diego Supercom- 
puter Center, typically using between 64 and 512 processors. Wall-clock time scales inversely 
with number of processors, indicating near-optimal parallelism; timing analyses are presented 
in Barranco & Marcus (2004). 

4. RESULTS OF NUMERICAL SIMULATIONS 

In the following presentation of the numerical simulations, we define a system of units 
such that Jl = H — p — 1, c s0 — QqH = 1, and RT = c 2 s0 = 1. With these 
units, the Mach number of a computed flow is e = max(-u) and the Rossby number is 
Ro = (1/2) max(cu 2 ). All of the simulations presented in this section were computed with 
T cooi —> oo. We will explore the effects of cooling in future work. 

4.1. Vortex initialization 

The key to initializing a 3D vortex is to guess a vorticity configuration that is as close to 
equilibrium as possible so that the initial accelerations are small and the flow can slowly relax 
to its true equilibrium; otherwise, the vortex will be ripped apart by large accelerations. As 
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discussed in Section 2.1, Moore & Saffman (1971) and Kida (1981) analytically determined 
exact steady-state solutions for 2D elliptical vortices of uniform vorticity embedded in a 
linear background shear. We use these 2D solutions as the starting point for constructing 
initial conditions for 3D vortices. We set v z and dv z /dt to be initially zero, which implies 
that the horizontal velocity fields at different heights are initially uncoupled and can be 
initialized independently. At each height, we construct a 2D elliptical vortex of constant 
Co z whose strength and shape satisfy equation (2-4) so that dv±/dt = 0. The ellipses at 
different heights do not necessarily have to have the same size, shape, or strength; in the 
following subsections, we will describe two different ways of "stacking" 2D elliptical vortices 
to create a 3D vortex. The anticyclonic core of the vortex is surrounded by a larger elliptical 
"halo" of weak cyclonic vorticity so that the net circulation in each horizontal plane equals 
zero: T = j A uo z dxdy = 0. We define a 2D streamfunction ip such that: vj_ = Vj_ x ipz and 
Cj z = —V 2 j_ip. Thus, once Q z is initialized at each height, we invert a 2D Poisson operator 
to obtain the streamfunction, which in turn generates the horizontal velocity field. The 2D 
pressure field can be found by taking the horizontal divergence of the horizontal momentum 
equation and inverting another 2D Poisson operator. The pressure will likely have a nonzero 
vertical gradient, which if left unbalanced, will lead to large vertical accelerations. We still 
have one more degree of freedom: we initialize the temperature field so that the buoyancy 
force exactly balances the vertical pressure force: 

f _ n_dm (4rl) 

Thus, all three components of the momentum equation are exactly balanced and the accel- 
erations are all initially zero: dv/dt = 0. Of course, only under very special circumstances 
would this procedure just so happen to also exactly balance the temperature equation. In 
general, dT/dt ^ initially, leading to an immediate evolution of the temperature field. 
Vertical motions will then be generated by the temperature changes through the buoyancy 
force, and these vertical velocities will then couple the horizontal motions at different heights. 



4.2. Tall, columnar vortex 

For our first initial condition, we constructed the simplest 3D analog of a 2D vortex: 
an "infinite" column created by stacking identical 2D elliptical vortices at every height. The 
vertical vorticity u z , horizontal velocity and enthalpy p/p were all independent of height, 
and the temperature perturbation was set identically to zero inside and outside the vortex. 
This initial condition was an exact equilibrium solution of the momentum and temperature 
equations: dv/dt = dT/dt = 0. However, it turned out that this quasi-2D vortex is unstable 
to small perturbations in 3D. 
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Ly = 3 




Fig. 3. — Tall, columnar vortex: y—z slices at x = of the ^-component of vorticity Q z . Blue 
= anticyclonic vorticity; red = cyclonic vorticity The initial elliptical column had aspect 
ratio x = 4 an d vorticity uj z = 0.625, corresponding to a Rossby number Ro = 0.3125. The 
times corresponding to each frame are: t/r or b = 0.0, 2.1, 4.2, 6.4, 51, 102, 153, 204. 
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Figure 3 shows the evolution of a perturbed columnar vortex. The domain dimensions 
for the simulation were (L x , L y , L z ) = (2,8,8), and the numbers of grid points/spectral 
coefficients along each direction were (N x , N y , N z ) = (64, 256, 256). The vertical velocity was 
forced to vanish on the top and bottom boundaries of the domain: v z (z = ±4) = 0. The 
minor and major axes (along the x and y directions, respectively) of the elliptical column 
were (A x ,Aj,) = (0.5,2.0), corresponding to an aspect ratio of x — 4. The magnitude of the 
anticyclonic vorticity in the core of the vortex was uj z = 0.625, corresponding to a Rossby 
number Ro = 0.3125. The maximum value of the initial horizontal speed was v± ;max = 0.12. 
In order to investigate the stability of the tall column, we did not stack the 2D ellipses exactly 
on top of one another, but instead offset the locations of their centers in the radial direction 
so that the centerline of the column had coordinates: {x c , y c } = {A sm(2iimz/L z ), 0}, where 
A = 0.1A r = 0.05 and m = 4. 

Because there is initially no vertical velocity or vertical enthalpy gradient, the different 
layers of the vortex do not communicate with each other, and tend to drift apart, carried 
by the ambient shear. As the vortex lines are stretched, a restoring force is generated which 
tries to vertically realign the vortex. The restoring force is partially successful in realigning 
a segment of the vortex around the midplane, but the parts of the vortex away from the 
midplane are sheared away, leaving a vertically truncated vortex straddling the midplane. 
Later, it appears that this truncated vortex goes through another instability and splits into 
two independent vortices above and below the midplane. These off-midplane vortices survive 
for well over a hundred orbits around the disk. 

From this simulation, it is clear that vertical stratification plays a key role in the vortex 
dynamics. The vertical component of the protostellar gravity is g z = Vt^z and the Brunt- 
Vaisala frequency is ujbv ~ ^oM/i?o- F &r from the midplane \z\ > H the gas in the disk is 
strongly stratified (Fr = tbv / T vor < 1), which inhibits the development of vertical motions 
which would couple the layers together. Thus, far from the midplane, the different layers of 
the vortex are free to drift apart, carried by the ambient shear. In contrast, there is virtually 
no stratification in the immediate vicinity of the midplane (Fr > 1), and a vertical flow 
rapidly develops to couple the layers together and prevent the ambient shear from ripping 
apart the vortex. 



4.3. Finite-height cylindrical vortex 

Motivated by the fact that tall, columnar vortices cannot maintain vertical coherence 
over more than a couple scale heights, we constructed a short, columnar vortex that was 
vertically centered on the disk midplane. We initialized a 2D ellipse in the midplane whose 
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ly = B 




Fig. 4. — Finite-height cylindrical vortex: y—z slices at x = of the z-component of vorticity 
Q z . Blue = anticyclonic vorticity, red = cyclonic vorticity The initial elliptical cylinder had 
aspect ratio x — 4 and vorticity uj z = 0.625, corresponding to a Rossby number Ro = 0.3125. 
The time between frames is At/r or .& ~ 60. 



Uc=i 




Fig. 5. — Finite-height cylindrical vortex: x—z slices at y = of the z-component of vorticity 
uj z . Blue = anticyclonic vorticity, red = cyclonic vorticity. The initial elliptical cylinder had 
aspect ratio \ = 4 and vorticity uj z = 0.625, corresponding to a Rossby number Ro = 0.3125. 
The time between frames is At/r or b ~ 60. 
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Fig. 6. — Finite-height cylindrical vortex: x—y slices at z = (first column) and z = 2 (second 
column) of the ^-component of vorticity u z . Blue = anticyclonic vorticity, red = cyclonic 
vorticity. The initial elliptical cylinder had aspect ratio x = 4 and vorticity uj z = 0.625, 
corresponding to a Rossby number Ro = 0.3125. The time between frames is At/r or & ~ 60. 
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Fig. 7. — Finite- height cylindrical vortex: y—z slices at x = of the temperature perturbation 
T. Blue = cool, red = warm. The initial elliptical cylinder had aspect ratio x = 4 and 
vorticity Q z = 0.625, corresponding to a Rossby number Ro = 0.3125. The time between 
frames is At/r orb m 60. 




Fig. 8. — Vortex lines for finite-height cylindrical vortex , at times t/r^ = 0, 170. 
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minor and major axes were (A x ,Aj,) = (0.5,2.0), corresponding to an aspect ratio of x = 4. 
The magnitude of the vorticity in the core of the vortex was Q z (z = 0) = 0.625 at the 
midplane, corresponding to a Rossby number Ro = 0.3125. The maximum value of the 
initial horizontal speed was v±^ max = 0.12. The vorticity was extended off the midplane 
according to a Gaussian profile: 

u z {z)=Co z {0)e~ z2,2H K (4-2) 

where the scale height of the vorticity was the same as the pressure scale height. For an 
anticyclone with Ro < 1, the inward Coriolis force is somewhat more dominant than than 
the outward centrifugal force, and the vortex must be a region of high pressure for horizontal 
equilibrium (see Figure 2). The high-pressure core extends only over a finite height, so that 
there is a vertical pressure force away from the midplane. In order for the vortex to be in 
vertical equilibrium, there must be cool, dense lids which provide a buoyancy force directed 
toward the midplane. 

Figures 4 and 5 show the time evolution of the ^-component of the vorticity in vertical 
slices y—z at x = and x—z at y — 0. Figure 6 shows the ^-component of the vorticity in two 
different horizontal slices x — y at z = (first column) and z = 2 (second column). Figure 7 
shows vertical slices y—z at x = of the temperature perturbation. The time between frames 
in all these figures is At/r or 5 m 60. These results were computed with the infinite vertical 
domain version of the simulation: the horizontal dimensions were (L x ,L y ) = (2,8), and the 
vertical mapping parameter was L z = 4. The numbers of spectral modes along each direction 
were (N x , N y , N z ) = (64,256,256). The three components of the momentum equation were 
exactly satisfied initially, but the energy equation was out of equilibrium from the start. The 
temperature field slowly evolved, generating a small vertical velocity which then coupled the 
horizontal layers together. After approximately a few dozen r or b, the vertically truncated 
vortex settled into a new, quasi-equilibrium (see frames 2-4 in Figures 4-7) that changed 
very little over the course of a couple hundred orbits through the disk. Figure 8 shows vortex 
lines for the initial condition and for the quasi-equilibrium steady-state at t/r or b = 170. 

We thought we had indeed found a stable steady-state, but surprisingly the vortex suf- 
fered a dramatic instability which ultimately resulted in the complete destruction of the 
vortex in the midplane (see frame 5 in Figures 4-7). The initial condition was symmetric 
with respect to the midplane, and the equations of motion (3-2) should have preserved this 
symmetry. However, the instability clearly had an antisymmetric component. We decom- 
posed the flow into its symmetric and antisymmetric parts. Figure 9 shows the maximum 
absolute value of the antisymmetric part of the ^-component of vorticity as a function of 
time. Initially, it was very close to zero, but grew from numerical round-off errors. Eventu- 
ally, a linear eigenmode emerged out of this numerical antisymmetric noise. The structure 
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of the mode preserved its spatial form for over ten orders of magnitude of growth, proving 
that it is indeed a linear instability. The e-folding time of the exponential growth was a few 

Torb ■ 

Thus, although the instability itself had a fast timescale, the vortex appeared to survive 
for a long time in its quasi-equilibrium state because a sufficient amount of numerical noise 
had to be generated in the antisymmetric modes to seed the eigenmode. In simulations in 
which we enforced symmetry about the midplane of the disk, the truncated vortex survived 
for longer periods of time, but eventually succumbed to slower symmetric instabilities. 

4.4. Off-midplane anticyclonic vortices and cyclonic azimuthal bands 

A closer inspection of frames 4-8 in Figures 4-7 revealed coherent regions of anticyclonic 
vorticity a couple scale heights above and below the midplane of the disk. These new 
vortices formed before the vortex in the midplane succumbed to the antisymmetric instability, 
and survived indefinitely in our simulations (over 400 orbits). The highly stratified disk 
supports neutrally stable internal gravity waves. As a midplane vortex slowly relaxes to 
its quasi-equilibrium, it oscillates and excites internal gravity waves which travel obliquely 
off the midplane of the disk. As the waves propagate off the midplane, conserving their 
flux of kinetic energy, the velocity and temperature perturbations increase, the waves break, 
and baroclinicity turns these perturbations into vorticity. The anticyclonic shear of the 
disk stretches the cyclonic vorticity into sheets or bands, whereas the anticyclonic vorticity 
perturbations roll-up into new, coherent vortices. Like the midplane-centered vortices, these 
off-midplane vortices are also high-pressure anticyclones. For a vortex in the upper half of 
the disk, the high-pressure core is balanced by a cool, dense top lid which pushes downward, 
and a warm, low-density bottom lid which pushes upward. These lids can easily be seen in 
frames 6-8 in Figure 7. The density and temperature of these lids are maintained by the 
vertical motions within the vortex: rising (falling) motion adiabatically cools (heats) the 
stably-stratified fluid. 

There are two intriguing features of these off-midplane vortices. First, the three vortices 
approach one another and strongly interact, stretching and squeezing each other. Clearly, 
these are robust features that are able to survive for hundreds of orbits despite strong per- 
turbations. One can see thin bands of cyclonic vorticity that keep the vortices separated; 
the bands are most prominent at closest approach (see the sixth frame down in column 2 
of Figure 6). Second, the off-midplane vortices appear to be "hollow" in the sense that the 
maximum vorticity does not occur at the center of the vortex, but on the rim. The Great Red 
Spot is probably the most famous example of a hollow vortex; the center of the GRS may, 
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in fact, be weakly counter-rotating (Marcus 1993). It is still a mystery as to how the GRS 
maintains its hollowness. Numerical simulations of 2D hollow vortices show that they are 
violently unstable; the low-vorticity fluid in the core switches places with the high-vorticity 
fluid on the rim, resulting in a vortex whose vorticity profile is centrally peaked (Youssef 
2000; Shetty et al. 2004). Here in our 3D simulations, hollow vortices are apparently stable. 

In order to get a better understanding on the structure of these off-midplane vortices, 
we isolated one of them and used it as a new initial condition. In this simulation, there is 
no vortex in the midplane, and except for the lone off-midplane vortex at z — 2, the rest of 
the domain is pure Keplerian with no initial pressure, temperature, or density perturbations 
anywhere else. Figure 10 shows the evolution of this isolated off-midplane vortex. In frame 

2 42r orb ), one can clearly see anticyclonic bands forming. By frame 4 126r orb ), the 
bands have rolled up and coalesced into two new vortices. The final state once again has 

3 off-midplane vortices. Most importantly, this simulation shows that that the vorticity 
for these new vortices did not come from the redistribution of vorticity remnants from the 
destruction of a midplane-centered vortex, but out of some new instability in the stratified 
region of the disk. In a future paper, we will focus on the mechanisms for vorticity generation 
and the sources of energy for these off-midplane vortices. 



4.5. Angular momentum transport 

One of the prime motivations for interest in coherent vortices in protoplanetary disks is 
whether they can transport sufficient angular momentum outward so that mass can continue 
to accrete onto the growing protostar. Following Balbus & Hawley (1998), one can relate the 
correlations of the fluctuating perturbation velocities (whether from turbulence or coherent 
vortices) to the eddy-viscosity model: 

«= (vxV y )/c 2 s , (4-3) 

where we have defined a density- weighted average (v x v y ) = f pv x v y dxdydz / f pdxdydz. For 
an isolated vortex that is fore-aft symmetric (e.g., v x is antisymmetric across the a;- axis, 
symmetric across the y-axis and v y is symmetric across the x-axis, antisymmetric across the 
y-axis) , one can show that the velocity correlation defined above vanishes; these symmetries 
must be broken in order to have non-zero transport. However, if the disk is filled with 
vortices, vortex-vortex interactions and mergers can break the fore-aft symmetry, allowing 
for angular momentum transport. 

In Figure 11, we plot the velocity correlations for the long-term simulation of a finite- 
height midplane vortex settling into a quasi-equilibrium, going unstable, and new, off- 
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midplane vortices forming. In the first panel, we plot the instantaneous correlations, which 
are highly time-dependent. Note that the correlations are both positive (outward transport 
of angular momentum) and negative (inward transport). The strong positive and negative 
spikes in this raw data are due to the interactions of the off-midplane vortices when they are 
at closest approach (see frame 6 in column 2 of Figure 6). In the second panel, we present 
a moving average (rectangular window with width 10r or &) of the same data. The smoothed 
correlations are always positive (outward transport of angular momentum) and correspond 
to an average a xs 1CT 5 . 

One must be careful when extrapolating the results of angular momentum transport in 
a local patch of the disk to the global transport in a whole disk. Our simulations are periodic 
in the azimuthal direction y, implying that there is a periodic array of vortices around the 
protostar, with azimuthal separation L y . However, vortices that are at the same radius 
tend to merge, potentially reducing the effective rate of angular momentum transport. The 
simulation shown in Figure 10 shows that new vortices readily form in unoccupied radial 
bands. Thus, there is a competition between mergers, which reduce the number of vortices 
per radial band, and the formation mechanism that repopulates empty bands. In future 
simulations, we will extend the azimuthal domain to investigate multiple vortices per radial 
band and explore this competition between mergers and formation of new vortices. 



5. SUMMARY & FUTURE WORK 

This is the first ever calculation of long-lived, robust 3D vortices in protoplanetary disks. 
The fundamental fluid dynamics are governed by three key timescales (see Table 1): r vor , 
the orbital period of gas around the vortex core; r or b, the orbital period of the vortex around 
the protostar; and tbv, the Brunt- Vaisala period, or period of vertical oscillations (e.g., 
internal gravity waves) in a stably-stratified atmosphere. These three timescales are roughly 
of the same order because they are all determined by the protostellar gravity. However, 
the stratification in the disk has a strong spatial dependence, making the dynamics more 
complicated. The disk can roughly be divided into two regimes: a weakly stratified region 
in the immediate vicinity of the midplane (\z\ < H , Fr = tbv/tvot > 1), and a strongly 
stratified region away from the midplane (\z\ > Hq, Fr < 1). Our numerical simulations 
have confirmed that stratification is a crucial ingredient for long-lived coherent vortices. 

Earlier, 2D studies of vortices in disks (e.g., Adams & Watkins 1995; Bracco et al. 
1998; Godon & Livio 1999, 2000) could not address the role of stratification. One way of 
interpreting their models is that a 2D vortex is a cross-section (in the midplane) of a tall 
column that penetrates through all the layers of the disk. Our numerical simulations (see 
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Figure 3) show that such columns, although they are exact steady-state equilibrium solutions 
to the 3D momentum and temperature equations, are unstable to 3D perturbations because 
stratification inhibits vertical motions which are needed to couple the layers together and 
lock the vertical alignment. Another possible interpretation of 2D analyses is that they 
capture the dynamics of vortices that are confined to the midplane of the disk. We initialized 
finite- height, cylindrical vortices straddling the midplane of the disk and allowed them to 
evolve into quasi-equilibria. These vortices were high-pressure anticyclones in which the 
inward Coriolis force balanced the outward pressure force. In the vertical direction, the 
vortex needed cool, dense lids to exert a buoyant force toward the midplane to balance the 
pressure force away from the midplane (see Figure 2). Long-term simulations of such vortices 
revealed that they too are unstable equilibria. We isolated an antisymmetric (with respect 
to the midplane) eigenmode that grew with an e-folding time of only a few orbital periods 
around the protostar. These vortices only seemed to last a long time because our initial 
conditions were symmetric with respect to the midplane, and the antisymmetric eigenmode 
was seeded only with numerical round-off error. 

We confess that when we started this work, we too believed vortices in protoplanetary 
disks would be either quasi-2D "infinite" columns, or short truncated cylindrical vortices 
confined to the midplane of the disk. Surprisingly, our simulations revealed an unexpected 
third possibility: off-midplane vortices a few scale heights above and below the midplane. 
These vortices are also high pressure anticyclones. For a vortex in the upper half of the 
disk, the high-pressure core is balanced by a cool, dense top lid which pushes downward, 
and a warm, low-density bottom lid which pushes upward. The density and temperature of 
these lids are maintained by the vertical motions within the vortex: rising (falling) motion 
adiabatically cools (heats) a stably-stratified fluid (see Figure 2b). Simulations with no 
vertical gravity and/or those with barotropic equations of state will fail to find such vortices. 
Also, numerical calculations must have sufficient resolution to capture the thermal lids, which 
typically have thicknesses much smaller than the local pressure scale height. 

Unlike the tall, columnar vortices and the finite-height cylindrical midplane vortices 
which we put in the disk by hand, the off-midplane vortices formed naturally from per- 
turbations in the disk. For example, as a midplane vortex oscillated and relaxed toward 
its quasi-equilibrium {not the cataclysmic motions as it succumbed to the antisymmetric 
instability), it excited internal gravity waves which propagated away from the midplane, 
steepened, and broke, creating vorticity (a baroclinic effect). Regions of anticyclonic vortic- 
ity rolled-up and coalesced into new anticyclones, whereas regions of cyclonic vorticity were 
stretched into thin azimuthal bands and sheets by the ambient shear. In simulations where 
we had only one off-midplane vortex as an initial condition, new vortices reformed, filling 
every available radial band. These off-midplane vortices are robust, surviving indefinitely 
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in our simulations (over 400 orbits) despite countless close encounters and interactions with 
other vortices. 

The symmetries of an isolated vortex preclude significant angular momentum transport. 
On the other hand, if the disk is filled with vortices that interact, merge and reform, those 
symmetries are broken and angular momentum can be transported at moderate rates. Thus, 
there is a competition between mergers of vortices that would reduce the number of vortices 
in a given radial band, and the formation of new vortices in unoccupied radial bands. 

In order to understand how vortices are created and maintained, and whether they are 
isolated or fill the disk, it is necessary to determine the sources of energy and vorticity. Klahr 
& Bodenheimer (2003) demonstrated that a globally unstable radial entropy gradient in a 
protoplanetary disk excited Rossby waves which also broke into large-scale vortices. Because 
the length scale of the entropy gradient was of order r, the horizontal scales of the vortices 
were much larger than the thickness of the disk, and so the fluid velocities were supersonic. 
Acoustic waves and shocks rapidly drain the kinetic energy from such vortices. We plan to 
extend the radial domain in our simulations and relax the shearing box boundary conditions 
so that we may investigate radial gradients as well. 

Coherent vortices located above and below the midplane will significantly affect the way 
dust settles into the midplane. Two-dimensional studies have shown that vortices are very 
efficient at capturing and concentrating dust particles (Barge & Sommeria 1995; Tanga et al. 
1996; de la Fuente Marcos & Barge 2001). If the vortices are located off the midplane, will 
the dust grains be inhibited from settling into a thin, dense sublayer? If the vortices have a 
significant vertical velocity, will they be able to sweep up grains that have already settled into 
the midplane? We propose that the trapping of dust in off-midplane vortices is analogous to 
the formation of hail in the Earth's atmosphere; turbulent vertical velocities prevent grains 
from falling out of the vortices until they have grown to some critical mass. In future work, 
we plan to incorporate simple Lagrangian tracking of dust grains (de la Fuente Marcos et al. 
2002), as well as two- fluid models (Cuzzi et al. 1993; Johansen et al. 2004). 
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Fig. 9. — Growth of unstable antisymmetric linear eigenmode. Because the initial condition 
was symmetric with respect to the midplane, the instability could not be manifested until 
round-off error populated the antisymmetric modes. Eventually, the eigenmode emerges out 
of the noise, and grows exponentially fast, with an e-folding time of a few r orb . 
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Fig. 10. — Simulation of an off-midplane vortex: x — y slice at z = 2 of the z-component 
of vorticity uj z . Blue = anticyclonic vorticity, red = cyclonic vorticity. The time between 
frames is At/r or .f, = 42. A single vortex from the last frame of Figure 6 was isolated and 
used as the initial condition for this run. Surprisingly, the other vortices reformed. 
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t/lorb t/x orb 



Fig. 11. — Velocity correlations as a function of time for the full simulation of a finite-height 
midplane vortex settling into a quasi-equilibrium, going unstable (at t ~ 200r orJ) ) , then new, 
off-midplane vortices forming. The first panel shows the raw velocity correlations, whereas 
the second panel shows the same data smoothed with a window of 10r or fe. 



